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1.0  INTRODUCTION 


The  quest  von  of  whether  or  not  the  wave  aberration  function 
of  an  isoplanatic  imaginq  system  is  uniquely  retrievable 
tron  the  system’s  point  spread  function  (PSD  appears  to  have  been 
fust  studied  by  A.  Wait  hot  and  n.  o*  Nei  1 1 ( 1 '  ~ } .  Basing  their 
s  lives*  i  oat  ions  on  the  well-known  result  that  the  so-called 
•tenor  a  1 1  zed  pupil  function 


!•  i  ».  .  ) 


ikoW(.,,  ) 

>° 


(1.1) 


terms  a  Fourier  transform  pair  with  the  coherent  spread  function 

-i  (  ;tX  +  i^v) 


( x ,  y  5  - 


h 


d,  F  (.»,?) 


e 


(1.2) 


Walther  and  O’Neill  correctly  nr>ted  that  the  problem  of  deducing 

w (  . , .  ■  f  rom  the  PSF 


h(x,y)a.-(x,y!^  (1. 3) 

really  amounts  to  that  of  retrieving  the  phase  of  ^(x,y)  from 
its  modulus.  The  problen  of  determining  the  phase  of  a  complex 
function  from  the  functions'  modulus  is  usually  termed  a  phase 
retrieval  problem  and,  as  discussed  by  Walther  and  O'Neill,  does 
not  generally  admit  a  unique  solution  unless  auxiliary  information 
of  some  sort  is  available  regarding  the  complex  quantity  whose 
phase  is  to  be  determined  (in  this  case  the  coherent  spread 
function  v(x,y)>.  As  an  example  of  this  indeterminacy,  we 


1.  e.L. O' Neill  and  A. Walther,  "The  question  of  phase  in  image 
formation,"  Opt.  Acta  10,  33-40  (1963). 

2.  A.  Walther,  "The  question  of  phase  retrieval  in  optics,*,  Opt. 
Acta  10,  41-49  (1963) . 


that  for  .my  ttv.intn:  system  |»« •hkohm mg  a  i.ulMlly  symmetric 
{'■‘Vil  Mnrt  ion,  the  wave  alu-irat  ton  (amnions  to  {  . ,  k  )  and 
Viola  identical  j’SK’a,  Thus,  f<»r  svn'h  systems,  estimates  of  the 
wave  aberration  function  baaoc*.  <  n  knowledge  only  of  the  PSF  will 
1h>  imlt't  et  mmit  e  at  least  up  t*-  t  rans !  o  rural  ions  of  the  type  given, 

I \ h  V«*  , 


l  : 


in;  the  phase  <■'  <x,yi  from  its  modulus  h*/*'{x,y) 
do.  iu>*  t :» *  t  the  phase  *>?  a  (unction  ( 1  ,i‘.  (  <x,y)) 

!:,nct  i'iii  an<!  the  modulus  of  its  Fourier 


i  :  t  he  r  -  m.-nt  let  as  suppose  that  besides  t  he  I'SF  we  also 

have  iv.u  lable  the  pupi  1  !  «m‘i  ion  ((,.  )  <ir,  what  amounts  t  o  t  h*. 

sure  *hin  t  he  r:o*;u!  is  It.,  )  <>t  tin*  ueiie  ru  1  1  ,s  ■<}  pupil  ("unction 

M<-.*.ins<-  ‘\,y>  amt  F  t  ..  1  t*ir  1  1  out  tei  transform  pair,  the 

.  ,  p  »  ,  .  >*’  •  ,  *1,  t.  '■**'•  ’  •  k.  •  *•.  t.  .....  •  ...  ...  *  ......  .  »  » .  RUt.ll.  I  111.  1%^/^ 

;■* -e  ernes  that  . > ' 

,vi'!i  tin 

‘  s  ar  sf  05  r  t  ;.*'.,  C  t 

l;  A  \ 

where  curtain 

workers  hive  investigated  the  possibility  of  deducin'}  the  phase  of 
1  coh*  ten*  i  y  I  i  !  .•.minuted  ob.»*ct  from  intensity  measurements 

pel  f  :  I 

1  •  r< ■;  •  . 

17,  .  1  .1  V  l  *  :»:• 

rot  r lev  e 


(i.  I ?  happens  that  precisely  this  later 

(  t  4 ) 

*l>ler  a;  1  ses  i:<  elect  1  mi  mio  os copy  * 


■  h  t  hr  lr.i-ie  and  diffraction  planes  of  an  electron 
*  (tin  '  .isni'ct  xm  <h*r<*hbotu  and  Saxton  J  devised 
tsuvf  r.ulf.  * f> l •  -  l  ast  Fe.urier  Transforms  (FFT's)  for 
1  f  .,  .n't  ion  from  its  modulus  and  the  modulus  of  its 


iM-.-r  transform.  Although  no  one  has  (to  our  knowledge)  been 
i!>  1  e  to  furnish  an  acceptable  uniqueness  theorem  as  to  the 
.elutions  ’blamed  by  the  derchberg-saxton  algorithm,  we 


Misell,  "An  examination  of  iterative  method  for  that 
s'datien  of  the  phase  problem  in  optics  and  electron  optics," 
Phys.  r>6.  2229-2225  (1973). 

a.j.  brenth  et  a!., The  problem  of  phase  retrieval  in  light  and 
electron  microscopy  of  strong  optics,  "Opt.  Acta  22,  615-62B 

(197*1. 

5.  r.w.  Oerchberg  and  w.o. Saxton,  "Phase  determination  from  imago 
and  diffraction  plane  pictures  in  the  electron  microscope," 
opt. ik.  14,  275-283  U971).  R.K. Gerchberg  and  W.Q.  Saxton, 

"A  practical  algorithm  for  the  determination  of  phase  from 
image  and  diffraction  plane  pictures,"  Opt  ik.  35,237-246  (.1972* 


2 


have 


a  a  others 


iiutu  5  n  ion  st  udies  that  the  algorithm  can 


bo  I’xpoiM c;;  to  yield  o.  unique  phase  reconstruction  so  long  as  the 
moduli  of  the  function  and  its;  Font  ier  i.  rants  form  are  adquatcly 


waits  >  i  ft! . 


The  Oerchbcrq-S.ixt.on  algorithm  is"  not  the  only  computation 
technique  available  (a;  deducing  the  phase  of  a  complex  function 
from  its  modulus  ami  the  modulus  of  its  Fourier  transform.  An 
alternative  scheme  has  also  been  proposed  by  Gerehbo rq  and  Saxton 
while  Gonsu 1 ves  provided  an  algorithm  specifically  tailored 


t<‘  the  :>roblc 


deduc  i  na  the  aberrations  of  a  Ions  from  the 


lens’  pupil  function  anal  joint  spread  function. 


The  preble"'  of  deducing  a  lens'  wave  aberration  function 
from  its  point  spread  function  i s  intimately  connected  to  the 
problem  of  ost i mat in«  wavefront  distortions  of  adaptive  optical 
imaging  systems  from  imagery  of  a  star  object.  In  this  later 
appl  i  eat  ion ,  the  image  of  the  star  object  is  detected  by  a  square 
array  of  photodiode:;  to  yield  an  estimate  of  the  PSP .  The  PSF 
estimate  is  then  input  to  a  digital  processor  programmed  with 
the  phase  retrieval  aluoi  The  output  from  the  processor  is 

A 

an.  estimate  K  ( -  )  of  the  wave  aberration  function  W(.-i»«)  which 

is  then  vised  to  determine  the  appropriate  control  signals  required 
to  drive  the  adaptive  imaging  system  to  a  state  having  a  more 


f>.  It. A.  Gonsalves,  "Phase  retrieval  from  modulus  data,"  J.Opt 
Soc.Ain.  6 6  ,  96  1  -‘MAI  ( 1976)  . 


favorable  aberration  function.  This  procedure  is  illustrated 
graphically  in  Figure  1.1. 


rieurc  1.1  Graphic  Illustration  of  *he  Use  of  a  Phase  Retrieval 
A 1  aor i thn  in  an  Adaptive  Imaging  Pa  tom. 


The  phase  retrieval  algorithms  discussed  above  can  also  be 
applied  in  situations  where  the  system  PSF  is  no*-  directly 
observable  but  can,  nevertheless,  be  estimate  .off  imagery  pro¬ 
duced  by  the  system.  For  example,  if  we  re?  t  our  attention  to 
a  single  isoplanatie  patch  the  image  interv-.  profile  l(x,y)  and 
object  intensity  profile  G(x,y)  are  rei  cd  v\.  the  equation 


1  ( x ,  y ) 


// 


o  ( x ' 


y  -  y 1 ) 


(1.4) 


where  h(x,y)  is  the  sytem  PSF  corresponding  to  the  particular 
isoplanatie  patch  in  question.  if  the  object  profile  0(x,y)  is 
known  the  liquation  (1.4)  can  be  inverted  (eg.,  by  use  of  Fourier 
transforms)  to  yield  the  PSF  h( . )  which  can  be  directly  input  to  a 
phase  retrieval  algorithm.  A  block  diagram  illustrating  the 
procedure  is  presented  in  Figure  1.2. 


; 


Figure  1.2  Block  diagram  or  control  loop  of  an  adaptive  imaging 
system  vising  an  "inverse  filter"  to  estimate  the  PSF  from 
imagery  of  a  known  object .  The  estimated  5'SF  is  input  to  a 
phase  retrieval  algorithm  for  determination  of  the  system 
aberrations  for  subsequent  control  of  the  optics. 


In  cases  where  the  scene  being  imaged  is  either  completely 
f.  or  partially  unknown  one  can  still  employ  the  aberration 

k  estimation  algorithm  discussed  above  in  con -junction  with  toch- 

r 

niques  which  adaptively  estimate  the  scene  being  imaged.  The 
overall  process  may,  however,  require  long  execution  times  and, 
i  thus,  may  not  even  converge  in  eases  where  the  image  motion  is 

|  pronounced.  For  this  reason  FIKONIX  has  developed  an  algorithm 


which  allows  one  to  estimate  an  optical  system’s  aberration 
function  from  two  images  of  an  unknown  scene  recorded  at  different 


system  do  see  us  values.  Ilona',  by  iho  use  of  dual  inane  planes 
such  as  shown  in  Fiqnre  1  .  t  in  oon  lunet.  io:i  with  the 
jfot'cmonl  ionod  ,i!-;on  i  hm,  it  is  posy  i  b  1  o  to  assess  the  wave 
abet 'Ml  ion  fuoo't  ion  of  t  he  i man  inn  system.  It  is  i  mpm  t  ant  to 
St  loss  that  no  a  ;>t  Lori  i  •  onnai  ion  rooi.i  r.i  1  no  the  scheme  be  inf) 
i .uevsecl  need  bo  known  and  no  dithcrinq  ol‘  elements  within  the 
adaptive  optical  system  is  required  to  deduce  the  wave  Trent 
distortions.  dust  as  was  i  he  case  when  the  scene  be  inn  i maced 
is  known  ,  the  ont  i  > pnenss  m.  asu  i  omen  t -wave  aberration 
est.  ima  t.  ion-correct,  ion  is  completely  nreomp !  i  shed  in  open  loop 
fashion  asina  only  the  imaqcry  recorded  <»'.’,•>  «  he  two  i maqe  planoa 

showti  in  !’  i  •  i  are  l  .  '■ . 


me..,  t  v 


beam  splitter 

im.iqo  plane  #1 
(CCD*  array) 


o  i : e  plan e  # 2 
l  ;  CCi)  array) 


*  i-isire  1  .  1  dassv  ;  i ,  n  system  with  adaptive  •.eeondary  (c.q.  , 
de:orr-able  r- i  r  ro r  i  ana  dual  lmnne  pl  ane:',  eon:,  t  t.  i  tv:  of  Charqe 
C:  ••  st.'  1  ed  Dev  l  ce  a  r  rays  .'OCI>>  located  at.  different:  distances  front 
second a  ry ,  !  r.  the  ibovo  f  i«ure  isn«v»e  plane  «1  is  located  at. 

best  focus  while  trauc  plane  *2  is  dot  oeused  by  the  (known) 
amount  of  <1  f  r td  me  *1. 


This  document  eon  fit.  i  tut  on  the  Final  Report  on  Contract 
F 10602- 77-c-Ol Tf  titled  "Analytical  Studies  of  Phase  Estimation 
Techniques."  The  objective  of  this  proqram  was  to  evaluate  the 
application  of  the  optica]  system  wavefront,  aberration  estimation 
techniques  discussed  above  to  multi-element  wide  field  of  view 
optical  con f ’ on  rat  ions .  To  this  end  it  was  necessary  to: 
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(I)  pert  nnn  an  analytical  i  nvost  i  gat  i  on  into  the  theory  of  wave 
aberration  estimation:  ( 2 )  develop  wave  aberration  estimation 
algorithms  to  be  used  in  computer  simulation  studies  anti,  (3)  to 
perform  computer  simulations  against  known  and  unknown  scenes  in 
the  presence  of  varying  amounts  of  noise  and  wavefront  error. 


The  estimation  foehn  i  qe.es  under  investigation  can  he 
classified  as  shown  m  Diagraml .  1 .  In  this  diagram  " 0 . S .  phase 
retrieval  algorithm"’  stands  tor  the  "Gerchbevg-Saxton"  algorithm 
which,  uses  the  mu.lt  iple  Fourier  transform  technique  described  in 
l  he  Svi>|K''‘  ’  ‘cv  ‘h>s  -it'trect  ami  shown  schematically  in 


Di  an  ram  1. '  Alternative  .System  Philosophies 


i'i'iui  i;  1.4  i'lock  diagram  i  1  1  i'.Kt  rat  in-?  the  operation 
Gorehborq-Saxton  (G.S.)  Algorithm  for  cst  i matin;'}  t  he 
the  general  I zod  pupil  function  !•’  .  from  its  modulus 

p  i  K 


or  the 
phase  ol: 


and  tho  modulus  «h<x  ,  v  ,  !  )  of 

r\  *  r*' 


its  Fourier  transform  . 

n ,  m 


Tho  results  for  the  analytical  investigation  phase  of  the 
program  are  reported  in  Chapter  2.  This  phase  of  the  effort  in- 
eluded  investigations  into  the  effects  of  polychromat  ic 
radiation,  detector  noise,  image  motion,  image  alias  inn, 
aberration  magnitude  and  spatial  frequency  content:  of  scenes  on 
the  process  of  phase  retrieval.  The  results  of  this  investigation 
ind  ieate  that  : 


1.  The  spectral  radiation  bands  anticipated  in  the 
HAI.O  mission  do  net  pose  a  problem  for  the  phase 
r  e  t  r  i  e  v  a  l  a  1  no  r  i  t.  hm . 

2.  Tho  NF.FD's  of  CCD  art  ays  anticipated  for  use  in  the 
HAI.O  mission  are  such  that  electronic  detector  noise 
will  not  seriously  degrade  tho  performance  of  the 
phase  retrieval  algorithms. 
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5.  Imago  aliasing  will  seriously  degrade  the 

performance  of  the  algorithms  but  can  be  overcome 
by  use  of  <i  multiplex  sump  l  im;  scheme  d  i  sous  sod 
in  Chapter  2. 

•1.  The  performance  of  the  phase  retrieval  algorithms 
is  optimum  for  moderate  amounts  of  aberration.  The 
performance  for  very  weak  and  very  strong  amounts  of 
a  be  r  rat  ion  de  pe  'ids  on  t  ho  s  i  q  n  a  1  to  n  o  i  s  e  r  a  t.  i  o 
and  on  the  spatial  frequency  content  of  the  scenes 
beino  i maced. 

5.  The  effects  of  scene  spatial  frequency  content  and 
finite  detector  size  on  the  phase  estimation 
algorithms  depend  op.  the  amount  of  aberration 
simulation  studies. 

A  large  part  of  the  work  performed  in  the  program  was,  of 
course,  devoted  to  actually  venerating  new  software  to  implement, 
the  various  abort  at  i«>n  ost.  imut  ion  techniques.  Software 
none  rated  during  the  course  of  the  program  includes: 


A  ;o:  -  o  T  t  •'  generate  tw< *-u i mens iona  i  Pel  s. 

A  program  to  perform  the  two-dimensional 
(*»«  rohber-.:-  Saxt.<.*n  al  aovithm. 

A  pros  ram  to  perform  the  ih>  v  a  n  e  y  - 1  ion  s  a  1  v  e  a  algorithm 
on  two-a  imens  tonal  data  from  the  one  imago  plane. 

A  s.roa  ran  to  per  form  t  ho  Oo  vaney -c.on  so  1  ves  alqorithm 
op.  two- :h  mens  vop.a  i  data  from  the  two  imago  plane. 

A  pro. vam  t.o  -venerate  polychrome*  1  c  one- 
u  i  men  s  i  on  a  !  OTP  •  y  , 


The  last  task  in  too  f-tatepont  of  Work  (Task  4.2.1)  calls 
computer  testinq  of  the  various  aberration  estimation 
ware  packages  developed  in  the  program.  The  results  from 
esentat i vc  simulations  are  presented  in  Chapters  .1,  4,  and  5. 


The  results  from  the  simulation  studies  using  the  Gerchbcrg 
Paxton  (G.S.)  algorithm  arc  presented  in  Chapter  3 .  This 
alcjonthm  did  not  perform  well  in  the  simulations.  In  particular, 
the  results  of  the  simulations  indicate  that  the  alqorithm  is  not 
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"robust"  in  two  dimensions  (i.e.,  its  performance  is  highly 
dependent  on  thr  choice  of  a  starting  cst. imat e)  .  because  ol  this 
it  was  decided  to  proceed  directly  to  the  two  dimensional 
Dovaney-Gon  sal  ves  algorithm  early  in  the  course  of  the  program. 

The  results  from  testing  of  the  two-dimensional  Devaney- 
(»on  salvos  algorithms  are  presented  in  Chapter  4.  We  say 
algorithms  because,  as  required  by  Task  4.1.6  of  the  Statement 
of  Work,  a  number  of  algorithms  were  developed  and  tested 
corresponding  to  different  models  of  the  wave  aberration  function 
the  principle  conclusion  to  be  drawn  from  the  simulation  studies 
presented  in  Chapter  4  is  that  these  algorithms  work  extremely 
well  and  are  very  robust;  so  long  as  there  is  no  limit.  t;o  the 
amount  of  computation  time  allowed.  Sometimes  extreme ly  long 
execution  times  are  required  by  the  algorithms.  The  reasons  for 
this  are  discussed  in  detail  in  Appendix  A,  along  with  a  state¬ 
ment  of  the  a  Igor  i  thin  and  .■«  more  complete  presentation  of  t:ho 
results  of  Chapter  4. 

Chanter  r>  presents  results  of  the  two  image  plane  simula¬ 
tions.  One  particular  wavefront  aberration  rs  used,  with  varying 
noise  levels,  detector  s.zes,  and  object  spectra.  Results  show 
that  the  algorithm  will  yield  good  estimates  even  in  the  presence 
of  largo  amounts  of  noise,  large  detectors,  and  large  objects. 
However,  the  same  problem  of  long  execution  times  is  present  as 
in  the  one  image  case.  Appendix  B  presents  a  ptopriot  rry  version 
of  Chapter  v> . 


0  AKAl.Vr.1S 


.'.1  SUMMARY 


In  this  Chapter  wo  review  t  ho  work  performed  on  the  tasks 
contained  with  the  “Analysis"  paragraph  (4.1)  of  the  “Statement 
of  Work"  for  she  Contract.  In  Section  2,2  we  discuss  the  effect 
ot  imdei  samp  1  i  n<j  of  an  image  of  a  star  object.  For  this  ease  it 
n>  shown  that  the  phase  retrieval  is  not  unique.  In  Section 
2. I  we  show  that  the  phase  retrieval  will  also  not  be  unique 
if  the  aberrations  arc  such  that  the  Optical  Transfer  Function 
(OTF)  is  well  approximated  by  the  so-called  "Geometrical  Optics" 
OTF .  Section  2.4  discusses  a  possible  technique  for  achieving 
Nyquist  rate  samp 1  mu  by  means  of  a  sampling  mult  iplex  procedure. 
The  effect,  of  finite  detector  size  and  of  scene  spatial  frequency 
content  are  discussed  in  Section  2.S  where  a  theorem  is 
established  which  shows  that  those  two  factors  have  precisely  the 
same  effect  on  phase  retrieval.  Finally,  Scot  1  on  2.  is  addresses 
the  question  of  the  effect  of  non -monochromatic  light  on  phase 
ret  rival . 


2.2  Tin:  F.FFFCT  OF  Aid  AS  l  NG  ON  PF.  ASF  RFTUIFVAl. 


For  the  case  of  a  star  object  the  phase  retrieval  problem 
reduces  to  that  of  determining  the  phase  of  the  coherent  spread 


iunct ion , 


v  ix) 


F  <«. ) 


-  i««x 


(2.1) 


from  its  modulus  J<,'(x)|  and  the  modulus  of  its  Fourier  transform 
F ( a i .  m  the  above  equation  «=  ^  u  is  a  "reduced  coordinate," 
with  f  being  the  focal  lonqth  of  the  imaging  system,  \  the 
wavelength  of  the  light  field  and  u  the  position  coordinate  over 
the  exit  pupil.  The  argument  of  tMx)  denotes  position  along 
the  i ma<>e  plane.  For  simplicity  of  presentation  we  have  restricted 


our  attention  t  o  t  he  oiio-u  intent*  i :>na  l  case.  The  results 
establ  i  Shod  below  hold  equally  well  in  the  case  of  two 
<1  i  men  <;  i  on ,  howovo  >  . 


We  note  that  since  the  integral  in  Equation  (2.1)  extends 
on  l  y  over  the  finite  interval  [-«,«]  wo  may  expand  the 
Generalised  Pupil  Function  F(<»)  in  a  Fourier  series  with 
period  2No  . 


!■■(.»)  X  C, 


i’Tw 


V-« 


(2.2) 


where  N  is  any  integer  '£  unity.  When  Equation  (2.2)  is  submitted 
into  Equation  (2.3)  one  obtains 


/ (» 0  /  ^ 

t 

<*>  *  E  c  I 

f  ,  ,  * (>>;r-  ' 

j  d  a  C  0 

V--  -  ) 

-!1 

0 

(2.  3) 

X 


V 


however,  from  Equation  (2.2)  we  conclude  that. 


C 

v 


! 

2fv 


o 


(2.4) 


where  we  have  used  Equation  (2.1).  Finally,  upon  substituting 
from  Equation  (2.4)  into  Equation  (2.3)  wo  obtain  the  saffpl inq 
series  representation  of  v(x): 
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It  is  important  to  note  that  tor  the  sampling  So  vies  (?,">) 
to  hold  the  quantity  N  must  bo  a  positive  integer  greater  than 
or  equal  to  unity.  It  N  =  1,  the  complex  Junction  T(x)  is  saic 
to  be  sampled  at  its  Nyqnist  Kate  while  for  N>1  it  is  sampled 
above  its  hyquist  Rate  (over  samp  loci)  .  Mow  the  modulus  square 
I V I  *  of  *i‘{x)  also  admits  a  sampling  series  representation.  In 
particular,  it  is  r.ot  difficult  to  show  that 


It  U) 


'  -  >  E  l  • 


M 


V  -  -  ... 


Ml 


(2  0i)l  *‘nC  o(X"  M(2o0)'-) 


(2.6) 


where  M  is  a  1  so  a  pos .  t.  » ve  nil  cnor  ureater  or  equal  to  unity. 

,  .  .? 

Again,  it  M  |  •;  ie|  *  s  said  to  be  sampled  at  its  Nyquist 

i  ate  while  i  !  M  -M  it  is  said  to  bo  oversamplod. 


Wo  note  that  the  Rygutst  rat.-*  tor  the  modulus  of  ^  (or, 

) 

equivalent  !y  !  »,  ( x  >  •  *' )  is  t-.-.ee  the  Ny-.uist  rate  for  v.  (This 
can  l.io  shown  !  o  be  cut'  to  the  t  act  .hat.  the  bandwidt  h  of  the 
modulus  of  t  t  .:nrt  ion  is  twice  ,  hat  of  the  function  itself.)  1 
particular,  the  maximum  sample  spacing  allowed  for  adequate 
sampling  of  Hx)  is 


\  1  ■ 


(2.7) 


whilt*  the  maximum  sample  spacing  allowed  for  sampling  of  the 


modulus  of  (x)  (o:  ,  equivalently,  of  .  (x)  “  1  is 


We  shall  now  show  that  if  the  modulus  of  ♦  is  sampled  at  or 
below  the  Nyquist  rate  for  t(i.e.,  if  Aj^|>  A^)  #  then  the 
retrieval  of  the  phase  of  n>from  its  modulus  and  the  modulus 
)F(a)l  of  its  Fourier  transform  will  not  be  unique  if  the 
pupil  function  is  a  constant. 


Theorem  I 

Let  the  pupil  function  f(a)  *|F(a)|  »  f^  be  constant. 
Further,  let  W(u)  be  a  continuous  but  otherwise  arbitrary 
wave  aberration  function  defined  over  the  exit  pupil  <_ 

-a  <_  a0  and  lot  w  be  a  real  number  in  this  same  interval.  Then 
any  aberration  function  of  the  form 

♦  Uq)  if  — ;iq  <  «<;  u 


Wu («)  * 


j 

|w  (a  -  M  -  (i 


(2.9) 


,)  if  u‘:  a  ■'  »  n 


will  generate  a  coherent  spread  function  (x)  havinq  the  same 

modulus  !  y  ( x )  |  at.  the  sample  points  x  *  n  A,  ;  n  »  0,  *1,  +2,. 

* 

Proof  of  Theorem  I 


From  Equation  (2.1)  we  have  that, 


lj»  (X)  a 


r°  /o 

f  F„l  "  c  '  ’of 

fy0 


, 

ir-w  (a)  -  v 

e  r  c  Au  .  U.10) 
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On  subst i tutiiu}  from  liquation  (2.9)  f  we  obtain 


wh  i  ch 


1  f  wo 

then  i 


*M  U) 


'•/ 


■>  i 


da  e 


i— -  W(u  -  n  ♦  <tQ) 


-  la  x 


0 

*  fo  I  <*« 


*T  wu 


•M  "  a0l 


-i«x 


(2.11) 


s  i  mpl  i  f  ios  t  o  booonc 


(x)  «  f0e'l(‘,’<'0)x 


/"- 


i  W  (a  *  ) 


-iaa* x 


-1! 


♦  f0o"l(ll4o0  )X|*  da’  o 


i j~  K(a') 


o 


-ici  *  x 


(2.12) 


put 


X  B  n  ‘  - 


•*0 


(2.11) 


0 


x  «  Jv  and  Equation  (2. 1 2)  (jives 


*0  i  «<a»)  (nA  ) 

V; 

(x~n  a  )  -  o  Ida*  o  o 


i  .  r 

7‘ 


e 


0 

(x  -n.  )  , 

V 


(2.14) 


(2.14) 


establishes  tho  theorem, 


l‘» 


which 


The  above  Theorem  shows  that  it  is  extremely  important  to 
sample  at  the  Nyquist  rate  appropriate  for  the  intensity 
of  the  image  in  order  to  obtain  a  unique  determination  of  the 
wave  aberration  function.  As  we  shall  show  in  the  following 
section  nonuniqueness  of  the  phase  retrieval  problem  is  also 
encountered  when  the  aberrations  of  the  system  are  so  severe 
that  geometrical  optics  gives  a  complete  description  of  the 
image  forming  properties  of  the  system. 


2.3  THE  EFFECT  OF  LARGE  ABERRATIONS  OF  PHASE  RETRIEVAL 


The  distinguishing  characteristic  of  a  highly  aberrated 
imaging  system  is  the  relative  unimportance  of  diffraction  and 
interference  effects  in  determining  its  image  forming  properties. 
In  particular,  it  is  common  practice  to  c‘  aracterize  such  systems 
by  the  so-called  "Geometrical  Optics  Optical  Transfer  Function" 
rather  than  the  exact  (within  the  limits  of  scalar  diffraction 
theory)  transfer  function  given  by 


x  y 


if 


OTF  (F  ,F  )  ■=  I  Idu  dv  f  (u,v)  f  (U+2\F*  F  ,  v42\Fi  F  » 


i  ~^|W(u,v)  -  W (u+2AF  ,  V42XF#  F  ) 

a  i*  ^  y  j 


(2.15) 


where  F#  is  the  F  number  of  the  imaging  system.  The  Geometrical 
Optics  Transfer  Function  is  obtained  from  the  exact  expression 
given  in  Equation  (2.15)  by  taking  the  Geometrical  Optics  limit 
•\  »  0.  The  details  of  performing  the  required  limiting 
operation  are  quite  straightforward,  (cf.,  K.  Miyamoto,  "Wave 
Optics  and  Geometrical  Optics  in  Optical  Design",  in  Progress 
in  Optics,  ed. ,  E.  Wolf  (North  Holland,  Amsterdam,  1966),  Vol.  I, 


p.43),  and  one  obtains 


OTFc..o.(IV  F 


>'  •// 


-i4irFi 


du  dv  c 
16 


f«,  3W(u,v)lt,  ?>W(u,v)l 
[x“u  *Fy - I 


(2.16) 


In  liquation  (2.16)  A  represents  the  area  of  the  unit  circle 
and  OTt'  is  the  Geometrical  Optics  OTF. 

In  the  following  analysis  we  employ  the  Geometric  Optics 
Transfer  Function  given  in  Equation  (2.16).  It  is  found  that 
this  quantity  by  it self  does  not  uniquely  determine  the  slopes 
<3 W/*u,  JW/9v)  of  the  wave  aberration  function.  In  other  words, 
within  the  framework  set  by  Geometrical  Optics,  wavefront 
estimation  via  the  use  of  phase  retrieval  techniques  is  not 
unique . 

It  is  extremely  important  to  note  that  this  nonuniqueness 
holds  only  in  the  geometric  optics  limit  when  diffraction  and 
interference  effects  can  be  completely  ignored.  These  effects 
will  be  present  no  matter  how  large  the  aberrations  and  it  is, 
evidently,  precisely  these  effects  that  are  responsible  for  the 
success  (and  apparent,  uniqueness)  of  phase  retrieval  techniques 
in  the  case  of  small  and  moderate  amounts  of  aberration.  If 
the  measurement  process  used  to  deduce  the  DTK  wore  perfect , 
then  the  effects  of  diffraction  and  interference  on  the  OTF  would 
be  observable  and,  hence,  the  process  of  phase  retrieval  would 
(presumably)  be  unique.  However,  when  the  aberrations  are  large 
and  measurement  noise  is  present,  these  effects  can  become 
"lost"  in  the  noise  with  the  result  that  the  Geometrical  Optics 
Transfer  Function  fits  the  noisy  data  quite  well  and  phase 
retrieval  will  not  lx?  unique. 

The  severity  of  the  aberrations  and  the  amount  of  measure¬ 
ment  noise  required  to  cause  phase  retrieval  techniques  to  fail 
is  not  known  at  present  and  can  probably  not  be  determined  in 
advance  even  for  the  simplest  optical  configurations.  However, 
preliminary  onc-dimensional  computer  simulations  indicate  that 
diffraction  effects  play  an  important  role  even  for  cases  of 
very  severe  (in  excess  of  ten  waves  of  aberration)  wavefront 
errors. 
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Finally,  it  should  oe  noted  that  although  the  results 
presented  below  establish  the  fact  that  the  wave  aberration 
function  cannot  be  uniquely  deduced  from  the  Geometrical 
Optics  Transfer  Function  alone  it  is  quite  conceivable  that 
certain  types  of  auxiliary  information  concerning  he  aberration 
function  might  lead  to  a  unique  specification  of  this  quantity. 
In  particular,  it  appears  lifcely  that  the  phase  retrieval  will 
yield  a  unique  solution  if  one  knows  in  advance  that  the  wave 
aberration  function  admits  a  finite  order  polynomial 
representation;  e.g.#  an  expansion  into  the  classical  wave 
aberration  polynomials. 

Proof  of  Nonuniqueness  in  the  Geometric  Optics  Limit 

For  simplicity,  let  us  consider  first  the  one-dimensional 
case  where  OTF^  *  reduces  to 

u  »U* 


with 


i  Fxg(u) 


OTFC.O.  <rx> 


du  c 


(2.17) 


=  -4n  F# 


3*»'(u) 


We  now  establish  the  following  Theorem. 


Theorem  I I 

If  g(u)  is  a  bounded  solution  to  the  integral  equation 
(2.17)  then  so  too  is  the  real  valued  bounded  function  d(u) 
which  satisfies  the  infinite  set  of  conditions 


/du  gn(u)  *  I 

•  ^  i 


du  d  (u)  ;  n  =  1,  2,  3,... 


(2.18) 


but  is  otherwise  arbitrary. 
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Proof  of  Theorem  II 


Since  the  integral  in  Equation  (2.17)  extends  over  a  finite 
interval  and  the  Taylor  series  expansion  of  exp[i  F  g(u)] 
is  uniformly  convergent  within  this  interval  (due  to  the 
assumed  boundedness  of  g(u))  we  may  write 


OTF 


G.O, 


<Fx>  - 


£  /‘ 

n=0  -1 


du  gn  (u) 


(2.19) 


If  the  bounded  function  d(u)  satisfies  the  infinite  set  of 
conditions  given  in  Equation  (2.18),  then  we  can  replace  each 
integral  on  the  r.h.s.  of  Equation  (2.19)  with  the  r.h.s. 
of  Equation  (2.1 6)  and  thus  obtein 


n  X 

OTPC.O.«V  *  £  —ITT-  f  '  *>  0"  ,u>  -/  dU  ' 

n-0  •'-1 


i  F  d (u) 
x 


(2.20) 


where  the  last  equality  results  from  the  assumed  boundedness  of 
d(u)  and,  hence,  the  uniform  convergence  of  the  sum.  This 
proves  the  theorem. 


Theorem  II  implies  that  the  quantity  g(u)  and,  hence, 
the  derivative  of  the  wave  aberration  function  jW(u)/3u  (i.e., 
the  wavefront  slope)  cannot  be  uniquely  deduced  from  the 
Geometric  Optics  Transfer  Function.  In  particular,  it  is  easily 
verified  that  the  infinite  set  of  conditions  given  in 
Equation  (2.18)  hold  if  segments  of  the  g(u)  function  are  inter- 
changed.  Even  if  one  imposes  the  further  condition  that  the 
function  <l(u)  be  continuous,  it  xb  still  not  difficult  to 
construct  a  function  which  satisfies  these  conditions. 

Figure  2.1  illustrates  one  method  of  generating  an  equivalent 
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function  d(u);  the  level  dQ  can  be  shifted  at  will,  showing 
that  there  are  an  infinite  number  of  such  functions  which 
yield  the  same  transfer  function,  .'rom  the  method  of 
construction ,  these  equivalent  functions  may  not  have 
continuous  fi*.st  derivatives. 


Figure  2.1  Two  Functions  Yielding  the  Same 
Geometric  Optics  Transfer  Function 
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An  equivalent  function  possessing  all  the  quantities  of  the 
original  function  g(u)  consists  of  the  “mirror  image"  of  the 
original  function.  That  is,  the  function  S(u)  =  g(»u)  as 
shown  in  Figure  2.2.,  will  have  the  same  transfer  function.  This 
can  be  proved  by  direct  substitution: 


Figure  2.2  Mirror  image  Functions 
Possessing  Identical  Transfer  Functions 


In  a  practical  situation,  the  wave  aberration  function  and 
hence  its  slope  g(u)  is  a  sampled,  discrete  process.  We  see 
that  if  there  are  N  independent  samples,  we  may  arrange  them  in 
any  order,  thus  generating  N!  different  functions  which 
correspond  to  the  same  transfer  function.  There  may  be  several 
sets  of  samples  which  assume  the  same  quantized  values;  thus 
the  number  of  distinguishable  displacement  functions  is 
accordingly  decreased,  but  this  is  not  important  to  this 
discussion.  The  important  point  to  note  is  the  lack  of 
uniqueness  in  determining  an  aberration  function  from  a 
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Geometrical  Optics  Transfer  Function. 


An  entirely  parallel  discussion  as  that  presented  above 
implies  for  the  two-dimensional  case.  In  particular,  we  have 
the  following  theorem; 

Theorem  III 


Let  3  W(u, v)/  0  u  and  3  W(u,v) /  3  v  be  the  first  partial 
derivatives  of  a  wave  aberration  function  which  yields  the 
Geometric  Optics  Transfer  Function. 


oftg.o.<fx'V 


-i4*Ft 

du  dv  e 


^.awiu.v) je.  3w(u,v)| 
F”Tu—  +Fy  5v  J 


(2.21) 


Further,  let  3W/?u  and  3W/3v  be  everywhere  bounded  within  the 
pupil  area  A.  Then  there  exists  at  least  one  other  wave 
aberration  function  W(u,v)  which  also  has  bounded  first  partial 
derivatives  3W/3u,  3W/3v  and  which  yields  the  same  Geometric 
Optics  Transfer  Function  as  does  W(u,v). 

We  omit  a  proof  of  Theorem  III  since  it  follows  directly 
from  the  discussion  presented  above  for  the  one-dimensional  case. 

2.4  A  PROCEDURE  FOR  INCREASING  THE  EFFECTIVE  SAMPLING  RATE  OF 

AN  ARRAY  OF  DETECTORS 

As  shown  in  Section  2.1  the  phase  retrieval  problem  does 
not  admit  a  unique  solution  if  the  PSF  is  sampled  below  the 
Nyquist  rate  of  (2/AFi)  samples  per  unit  length.  Throughout  this 
document  it  is  assumed  that  this  sampling  requirement  is  met. 
However,  it  is  anticipated  that  the  actual  detectors  to  be  used 
in  the  HALO  mission  will  be  too  large  to  meet  this  rather 
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stringent  sampling  requirement  and,  thus,  it  was  necessary  to 
explore  alternative  procedures  for  achieving  the  required 
sampling  rates. 

The  procedure  described  in  this  section  uses  a  relative 
image  motion  between  the  detector  array  and  the  image  from 
which  the  aberration  estimation  is  being  performed  to  obtain 
a  higher  sampling  rate  than  would  normally  be  achievable  with 
the  detector  array.  Since  the  effect  of  the  finite  areas  of 
the  detector  surfaces  simply  modifies  the  spatial  frequency 
spectrum  of  the  object  being  imaged  (see  Theorem  IV  in 
Section  2.4)  we  shall,  in  the  following  discussion,  limit  our 
attention  to  an  array  of  point  detectors.  The  results  obtained 
hold  equally  well  for  the  case  of  detectors  having  finite  areas. 

The  proposed  technique  is  most  easily  introduced  in  a 
one -dimensional  context..  In  Figure  2.3  is  shown  a  one- 
dimcnsional  PSF  which  is  sampled  by  the  linear  array  of  point 
detectors  having  a  sample  spacing  of  A  .  It  is  assumed  that  the 
PSF  moves  across  the  detector  array  with  a  uniform  velocity  V 
and  that  it  does  not  change  shape  over  the  time  interval  (A  /V) . 
At  time  tj,  the  detectors  will  measure  the  values  of  the  PSF 
labeled  by  the  number  one  in  the  Figure  (The  detectors  are 
assumed  to  be  "turned  on"  only  in  a  very  short  time  interval 
centered  at  tj).  At  time  ^2=tl  *  d/v~t^+  A /2v  the  detectors 
will  measure  the  values  of  the  PSF  labeled  by  the  number  two  in 
Figure  2.3.  Decause  the  time  increment  between  the  reading 
at  1 1  and  t2  was  selected  to  be  (A/2v)  the  above  process  is 
seen  to  yield  equally  spaced  samples  of  the  PSF  with  the 
samples  being  A /2  units  apart.  It  is  clear  that  by  choosing 
smaller  values  of  the  time  increment  even  hiqher  sampling  rates 
can  result. 
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Array  Axis 


Figure  2.3  Two  sets  of  samples  of  a  PSF  obtained  at  the  two 
times  t^  and  t2  =  ♦  d/v.  A  relative  linear  velocity  of  V 

is  assumed  to  exist  between  the  PSP  and  the  detector  array. 


An  entirely  analogous  situation  as  that  described  above 
results  in  the  two-dimensional  case.  In  two  dimensions,  however, 
it  is  important  that  the  motion  of  the  PSF  at  an  angle 
relative  to  one  of  the  detector  axes.  This  condition  is  necessary 
in  order  to  achieve  a  high  sampling  rate  in  two  orthogonal 
directions. 
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In  Figure  2.4  is  shown  an  array  of  point  detectors  (large 
dots)  having  a  sample  spacing  of  A.  Assume  that  a  PSF  were 
to  sweep  across  this  array  at  a  velocity  of  V  and  at  an 
angle  at  45°  with  respect  to  the  array's  horizontal  axis.  If 
the  detector  array  were  to  make  two  successive  instantaneous 
measurements  separated  by  a  time  interval  of  (a//5TV)  seconds 
apart  then  the  two  sets  of  measurement  together  would  yield 
samples  of  the  PSP  at  all  sample  points  shown  in  the  Figure. 
(The  large  dots  would  be  the  sample  positions  at  one  time  and 
the  smaller  dots  the  sample  positions  at  the  other  time.)  The 
overall  process  would  thus  be  equivalent  to  sampling  the 
stationary  PSF  by  a  regular  square  array  of  point  detectors 
having  the  geometry  3hown  in  the  box  in  Figure  2.4. 


Figure  2.4  Samples  of  a  rsr  moving  at  a  uniform  velocity 
at  45°  relative  to  the  horizontal  axis  of  an  array  of 
point  detectors  (large  dots). 
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A  higher  sampling  rate  than  that  shown  in  Figure  2.4  can 
be  achieved  if  the  angle  0  between  the  array  axis  and  the 
direction  of  motion  of  the  PSF  is  smaller  than  45°.  However, 
it  is  easily  shown  that  in  order  to  obtain  an  equally  spaced 
square  matrix  of  samples  this  angle  can  only  assume  values  which 
satisfy  the  equation 

Tan  o*  1/n;  n  «  1#  2,  3,  .  .  •  (2,22) 

A  plot  of  the  allowable  values  of  ©  versus  n  is  shown  in 
Figure  2.5. 
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Figure  2.5  Values  of  0  Satisfying 
Equation  (2.22) 


t~ . .  - T-  I  ■■■  —  .y— ■  «■  I  -■■■—> 

1  ^  d  t  t>  7 


26 


The  condition  ( I'q .  2.22)  results  from  the  requirement 
that  the  relative  motion  between  the  detector  array  and  the 
PSK  bo  such  that  at  the  beginning  of  a  new  sampling  sequence 
the  relative  position  of  the  PSF  and  the  detector  array  be  the 
same  as  it  was  at  the  initiation  of  the  previous  sampling 
sequence.  In  the  case  shown  in  Figure  2.4  the  sampling  sequence 
consisted  of  two  9ets  of  measurements.  At  the  initiation  of 
the  sequence  the  measurement  points  are  indicated-by  the  large 
dots  while  at  the  second  sampling  time  the  measurement  points 
would  be  the  small  dots.  Since  the  direction  of  motion  is 
45°  the  next  sampling  time  (assuming  samples  to  be  taken 
(i\\j  2v)  seconds  apart)  would  occur  again  at  the  locations  of 
the  large  dots  and,  thus,  the  whole  process  would  bo  repeated. 

In  Figure  2.6  is  shown  the  case  o it  n  =  2  or  (>=26,56°.  The 
sampling  sequence  consists  of  five  equally  spaced  measurements 

occurring  (1A/  5  A/v)  second  apart.  (in  general,  the  number 

”  2 

of  measurements  in  a  sampling  sequence  is  n  ♦  1  and  the  time 

between  samples  is  thus  ((l/ynZ  ♦  1)  A/v  seconds).  At  the 
sixth  measurement  time  the  relative  position  between  the  sampling 
array  and  the  PSF  is  the  same  as  at  the  initiation  of  the  sampling 
sequence  sc  the  whole  process  repeats  after  every  five  measure¬ 
ments. 


'-?r? 
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One  of  the  major  problems  encountered  will!  sampled  <lata 


imaging  systems  is  aliasing  introduced  by  undersampl ing  in  the 
image  plane.  In  addition,  the  finite  si'/e  ol  each  element  of 
an  array  of  detectors  has  the  effect  of  decreasing  the  object 
spectrum  at  high  spatial  frequencies  and,  thus,  leads  to  small 
signal  to  noise  ratios  at  high  spatial  frequencies,  both  of 
these  effects  have  an  impact  on  the  process  ol  aberration 
function  estimation.  In  tho  following  discussion  we  present, 
the  analysis  which  relates  the  output  of  a  two-dimensional 
array  of  detectors  to  the  image  being  detected  and  in  the  process 
establish  the  result  that  scene  spatial  frequency  content  and 
finite  detector  size  have  precisely  tho  same  impact  on  phase 
retrieval . 


We  shall  assume  that,  the  image  l(x,y)  of  an  object.  0 ( x ,  y ) 
is  sampled  by  a  regular  array  of  identical  detectors  such  as 
shown  in  Figure  2.7,  Tn  this  Figure 
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(2.23) 


where  l/2f  is  the  center  to  center  separation  of  the  detectors 

o 

and  n  and  ro  are  integers. 
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Figure  2.7  Sampling  Array  Geometry 
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in  the  absence  of  noise  the  response  of  each  detector  will  be 
assumed  to  be  given  by 


ln =//d*' 


dy'  I(x',  y’)  D(x*-xn,  y*  -  yR) 


(2.24) 


where  D(x',y')  *  0  if  jx’lor  {y'|>  1/4  f  . 


By  substituting  the  Fourier  integral  representation 


I  (x\  y* )  » 


i  {[  -  lfv,+v1 

“5  }]  d\dvy  1  <W  e  L  ,2.25, 


into  Equation  (2.24)  we  find  that  i  _  becomes 

n,m 


i  /Id K  dK 

n'B  <2n  )2JJ  *  y 


1  (Kx»Ky)  D(Kx,Ky)  e 


Kv  *ol 


(2,26) 


where 


dx'dy*  D ( x ' ,  y‘)  e 


-i[v  ,  v] 


(2.27) 


and 


■If 


I (K  «  K  )  »  f  |dx ’ dy •  lix'.y'le 
x  y 


-i  fv'  ♦  V'J 


=  0(KX.  Ky)  OTF  (Kx,  Ky) 


(2.28) 


with  0  (K  ,K  )  being  the  spatial  frequency  spectrum  of  the 
x  y 

object  being  imaged  and  OTF  (K  ,  K  )  the  Optical  Transfer 

x  y 

Function  of  the  imaging  system.  Finally,  on  substituting 
Equation  (2.28)  into  Equation  (2.26)  we  find  that 
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dKxdKy  o(Kx,  Ky)  OTF(Kx>  Ky)  e 


1  [Kx5? *  Ky*r?  1 
*-  o  o 


(2.29) 


where 


o(Kx,  Ky)  «  D(Kx,  Ky)  0(Kx,  Ky)  (2.10) 

We  conclude  from  Equation  (2.29)  that  the  effects  of  finite 

detector  size  as  manifested  in  the  detector  frequency  response 

function  D (K  ,  K  )  and  object  spatial  frequency  content 
v  y 

0(Kx#  K  )  have  entirely  equivalent  effects  on  the  output  in  m 
from  an  element  of  the  detector  array.  This  "equivalence 
property"  allows  us  to  conceptually  replace  an  actual  detector 
array  by  an  array  of  point  detectors  so  long  as  we  also 
replace  the  actual  object  spectrum  O ( ,  Ky)  by  the  modified 
spectrum  0(K  ,  K  _)  as  defined  in  EQuation  (2.30).  This 

*  x  y  ~  ■■  . . . 

equivalence  is  stated  precisely  in  the  following  Theorem: 

Theorem  IV 

I.et  the  image  of  an  object  having  a  spatial  frequency 
spectrum  0(K  ,  K r)  bo  sampled  by  a  regular  array  of  identical 

detectors  having  a  center  to  center  separation  of  l/2f  and 

\-  w 

frequency  response  funct  ion  l)(K  ,  K  )  .  Then  the  output  from  this 

x  y 

detector  array  is  precisely  the  same  as  would  be  obtained  from 

detecting  the  image  of  an  object  having  the  spatial  frequency 

spectrum  o(K  ,  K  =  0(K  ,  K  )  D(K  ,  K  )  by  a  regular  array  of 
x  y  x  y  x  y 

point  detectors  having  the  same  center  to  center  separation 
U.c.,  l/2fo  ). 


31 


Proof  of  Theorem 


The  proof  of  the  Theorem  follows  at  once  from  Equations 
(2.29)  and  (2.30). 

Comments 

The  above  result  is  extremely  important  in  that  it  allows 
us  to  always  deal  with  an  array  of  point  detectors  without  any 
loss  of  generality.  Because  of  this,  for  example,  the  analysis 
presented  in  the  proceeding  section  is  directly  applicable  to 
arrays  of  non-point  detectors.  In  addition,  it  shows  that  the 
problem  of  phase  retrieval  from  a  star  object  using  a  regular 
array  of  non-point  detectors  is  precisely  the  same  as  the 
problem  of  phase  retrieval  from  non-point  objects  using  an 
array  of  point  detectors.  This  later  equivalence  turns  out  to 
be  extremely  important  in  the  simulation  studies  presented  in 
Chapter  4. 

2.6  THE  EFFECT  OF  NON-MONOCHROMATIC  RADIATION  ON  THE  POINT 

SPREAD  FUNCTION 

In  order  to  determine  the  effect  of  poiychromat ic  radiation 
on  the  process  of  phase  retrieval,  a  number  of  (one-dimensional) 
multi-spectral  PSF's  and  OTF's  were  generated,  having  various 
amounts  of  coma  and  do focus.  The  input  spectra  for  these  PSF's 
and  OTF's  contained  three  to  five  wavelengths  and  had  a  total 
width  of  10%  to  20%  of  the  central  vavelenth.  The  results 
indicate  that  the  introduction  of  many  spectral  components 
affects  the  detailed  structure  of  the  PSF  but  not  its  general 
shape.  Another  way  of  putting  this  is  that  the  OTF  is  affected 
mainly  at  the  higher  spatial  frequencies. 
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Method  of  Calculation 

Let  P(x,  A)  be  the  point  spread  function  at  wavelength  at 
point  x.  Then  the  multi-spectral  PSF  is  given  by 


P(x) 


P(x, 


(2.31) 


where  the  quantities 
different  wavelengths. 


represent  the  relative  weights  of  the 
For  each  wavelength  P(x,a  )  is  given  by 


where  F#  =  the  F  number  of  the  system 

a  -  the  aperture  coordinate/the  aperture  radius 

f(a)  =  the  pupil  function  of  the  system 

W(f»2  ~  the  total  phase  shift  due  to  the  wave  aberration 
y  function  W (<» )  at  wavelength  A 

It  is  important  to  note  that  it  is  tactily  assumed  in  this 
analysis  that  the  wave  aberration  function  W(u)  is  independent 
of  wavelength.  Such  an  assumption  is  valid  for  optical 
imaging  systems  employing  high  quality  mirrors  such  as  is 
envisioned  for  HALO. 


The  plots  presented  in  Figure  2.8  through  2.16  show  the 
total  phase  shift  at  the  central  wavelength;  i.e., 

«£!  >  E  A  (2. 33 

A  n  n 

the  multi-spcctral  (polychromatic)  PSF,  and  the  modulus  and 
phase  of  the  polychromatic  OTF  for  various  choices  of  the 
expansion  coefficients  (aberrations)  A^  and  for  various 
wavelength  spectra.  Figures  2.8  -  2.10  show  line  spread 
functions  and  OTF* s  calculated  for  an  aberration  of  one  wave  of 
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Figure  2. 10  Polychromatic  Case  For  One  Wave  of  Defocus 


Single  Wavelength 


Figure  2. 13  Polychromatic  Case  For  One  Wave  of  Coma 


Polychromatic  Case  For  Oct  Wave  ol  Defocus  Plus  Oie  Wave  of  Coma 


I1  In  H  i' i: 


defocus.  The  first  figure  ahows  the  monochromatic  case, 

while  the  next  two  figures  show  the  case  of  three  wavelengths 

with  a  total  spread  oif  -*-10%,  and  five  wavelengths  with  a  total 

spread  of  +20%.  Figures  2.11  -  2.13  and  Figures  2.14  -2.16 

show  a  similar  series  for  one  wave  of  coma,  and  for  one  wave 

of  coma  plus  one  wave  of  do focus. 

*  ■ 

hs  mentioned  above,  these  results  indicate  that  the 
polychromatic  and  monochromatic  OTF’s  are  quite  similar, 
differing  mainly  at  high  spatial  f requoncies.  th  this  sense  the 
effect  of  polychromat ic  radiation  on  the  process  of  wave 
aberration  estimation  is  quite  similar  to  the  effect  of  detector 
noise  on  this  process.  Because  the  phase  retrieval  algorithms 
appear  to  be  highly  insensitive  to  moderate  amounts  of  detector 
noise,  it  is  reasonable  to  conclude  that  they  will  likewise  be 
insensitive  to  the  effect  of  moderate  amounts  of  spectral 
broadening  in  the  detected  radiation.  This  hypothesis  can, 
of  course,  be  tested  in  computer  simulations  although  it  was 
decided  that  such  simulations  were  not  warranted  within  the 
current  effort. 


3.0  THE  GERCHHERC-SAXTON  ALGORITHM 


The  Gerchberg-Saxton  (G.S.)  algorithm  is  a  method  which 

permits  determination  of  wavefront  aberrations  from  measured 

point  spread  functions  by  the  use  of  repeated  Fourier  transforms. 

The  overall  operation  of  the  algorithm  is  shown  in  Figure  1.4 

of  Chapter  1  and  in  more  detail  in  Figure  3.1  below.  It  basically 

consists  of  going  back  and  forth  between  the  generalized  pupil 

function  and  the  coherent  spread  function  by  means  of  Fourier 

transforms,  each  time  keeping  the  phase  part  of  the  function  and 

replacing  the  magnitude  by  the  known  correct  magnitude,  which 

is  the  square  root  of  the  measured  point  spread  function  for 

the  coherent  spread  function,  and  unity  for  the  generalized 

pupil  function.  There  is  no  proof  that  the  algorithm  must 

converge,  or  that  if  it  does,  it  will  converge  to  the  correct 

aberration  function.  Clearly  one  possible  mode  of  convergence, 

which  is  the  desired  result,  is  that  W  becomes  equal  to  the 

2 

actual  wave  aberration,  and  that  |CSF|  is  exactly  equal  to  the 
measured  h.  If  a  different  W  were  found  such  that  we  still  had 
|  CSP|  =  h,  then  this  would  be  an  example  of  non-uniqueness. 

This  has  never  occurred  in  our  simulations,  although  we  cannot 
prove  its  impossibility.  A  third  possibility  of  convergence 
is  one  in  which  neither  W  nor  CSF  over  reach  the  proper  values, 
but  the  estimates  remain  unchanged  each  time  around  the  loop. 

This  last  possibility  in  fact  seems  to  be  a  coimnon  occurrence 
although  we  are  unable  to  state  if  the  algorithm  has  truly  con¬ 
verged  to  a  poor  result,  or  is  merely  changing  at  a  very  slow 
rate . 

Previous  work  performed  prior  to  this  contract  extensively 
tested  the  G.S.  algorithm  in  the  one-dimensional  case  with 

(7) 

excellent,  results.  No  examples  of  false  convergence  were 

T~.  A .  J .  De vane j .  H .  A .  Gon s a  1  ve s  and  R.Chidlaw,  "Appl  ication  of  phase 
retrieval  techniques  to  adaptive  imaqino  systems,"  J.  Opt. 

Soc.  Am.  ( A) 6  7 ,  1422  (3977). 


Gerchberg-Saxton  Phase  Retrieval  Algorithm 


observed.  The  two-dimensional  G.S.  algorithm  was  tested  on 
this  contract  for  a  number  of  cases.  For  some  initial  (random) 
estimates  of  the  phase,  the  procedure  worked  well.  For  others, 
convergence  is  extremely  slow  (perhaps  on  the  order  of  10,000 
iterations)  if  indeed  it  would  even  be  reached.  A  convergence 

( g  \ 

acceleration  technique  inspired  by  Pienup  yielded  marginal 
improvement  in  both  the  fast  converging  and  slow  converging 
cases;  but  the  slow  converging  cases  remained  too  slow  for  use. 
The  fast  converging  cases  can  be  driver,  into  instability  by  too 
much  of  this  procedure. 

Table  3.1  is  a  list  of  results.  The  amount  of  aberration  is 

2 

expressed  in  radians,  in  terms  of:  defocus  -  P,r  ,  third-order 

4  J  2 

spherical  -  F^r  ,  and  third-order  astigmatism  -  P^r  cos  9  . 

The  initial  estimate  of  the  phase  aberration  is  chosen  as  an 

array  of  random  numbers.  The  RMS  value  of  the  array  is 

adjusted  depending  on  the  Strohl  ratio  of  the  PSP.  A  linear 

tilt  is  added  if  the  peak  value  of  the  PSF  is  not  at  the  array 

center.  The  random  numbers  are  generated  by  a  system  routine 

on  the  PDF  11/70  called  RAN.  It  requires  two  random  number 

"seeds"  to  determine  a  particular  pseudo-random  sequence.  Three 

sets  of  seeds  were  used,  as  listed  in  Table  3.2. 


o.  J.R.  rienup,  "Reconstruction  of  an  object  from  the  modulus 
of  its  Fourier  transform",  J.  Opt.  Soc.  Am.  (A)  67,  1389  (1977)/ 
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Aberration 

(radians) 

Random 

Number 

seeds 

Itera¬ 

tions 

Per¬ 

formed 

N ,  N2 

mu 

Results 

Start  Min.  Finish 

26 

■1 

X 

3.7  8.3 

27 

i  4 

B 

Q 

XII 

315.0  112.0  115.0 

28 

P3  *  2) 

A 

Wm 

16,32 

bad 

29 

P4  =  4 

100 

60  0.03 

bad 

30 

(r5  *  S) 

300 

1 

X 

402  to  133  (de¬ 
creasing  very 
slowly) 

31 

40 

XV 

402  to  178  (bad) 

32 

40 

k 

XIII 

(very 

strong) 

402  to  169  (bad) 
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TABLE  3.2  RANDOM  NUMBER  SEEDS 


A  11427  9223 

B  2689  8957 

C  4567  13879 

The  lens  aperture  is  a  circle  of  maximum  diameter  which 
will  fit  into  an  N  x  N  array  of  points.  This  is  buffered  out  to 
N2  x  N2  with  zeros.  If  N2  is  twice  N,  the  resulting  PSF  is 
sampled  at  the  Nyquist  rate. 

The  column  labeled  "acceleration”  refers  to  the  parameters 
using  in  attempting  to  accelerate  convergence.  If  left  blank, 
the  original  G.  S.  algorithm  was  used.  In  the  "results" 
column,  a  judgment  of  how  well  the  first  aberration  estimate 
compares  with  the  actual  input  aberration  is  listed.  Eventually 
a  measure  of  the  error  between  the  original  and  estimated  PSF's 
was  added  to  provide  a  more  exact  measure. 

Cases  1,  2,  and  3  show  the  effect  of  increasing  the  number 
of  iterations  of  the  G.S.  algorithm,  which  is  an  improvement 
in  the  estimate  of  the  phase  aberration.  A  different  random 
number  seed  yields  an  estimate  that  bears  no  relation  to  the 
actual  aberration,  in  Case  4.  Other  results  show  that  this  case 
has  not  really  converged  but  is  slowly  changing. 

Cases  12,  13,  and  14  have  a  smaller  amount  of  aberration, 
with  the  random  number  seed  in  Case  12,  after  only  16  iterations, 
the  final  estimate  is  definitely  headed  for  convergence  at  the 
actual  starting  aberration.  With  different  random  numbers, 
even  80  iterations  yields  bad  results.  In  Case  11,  the  effective 
sampling  of  the  point  spread  function  has  been  doubled,  by 
reducing  the  number  of  points  over  the  aperture  to  8,  but  this 
case  still  failed  to  converge. 
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Cases  5  th"ough  10  illustrate  the  first  attempt  at 
accelerating  the  algorithm.  The  G.  S.  algorithm  was  modified 
by  the  addition  of  a  term  F •  (  /h  -  CSF  ) ,  in  the  box  labelled 
"modified"  in  Figure  3.1.  This  term  causes  the  CSF  to  be 
corrected  not  to  the  known  CSF  modulus*  /h,  but  rather  shifted 
by  an  amount  proportional  to  the  difference  between  the 
estimate,  ICSF|  ,  and  the  desired  quantity  /fi,  but  in  the 
opposite  direction.  The  parameter  F  denotes  the  strength  of 
this  correction.  This  new  term  has  not  altered  what  may  be 
termed  an  eigenfunction  of  the  loop;  if  a  W  which  exactly 
yields  a  point  spread  function  h  is  put  into  the  loop,  it  will 
return  unchanged.  Thus  any  function  W  at  which  the  old  algorithm 
would  converge  will  still  be  such  a  function  in  the  new  algorithm. 
It  was  hoped  that  this  would  eliminate  the  false  convergence 
problems. 

In  Cases  3  through  8,  F  was  set  equal  to  the  value  specified 
in  the  table  from  the  beginning  of  the  algorithm.  For  low  F*s, 
this  converged;  for  large  F's,  the  algorithm  did  not.  By 
waiting  for  10  iterations  before*  applying  an  F  of  0.05,  there 
is  still  no  convergence  in  Case  9.  But  in  Case  10,  applying 
F  =  0.05  after  15  iterations,  the  algorithm  does  converge.  It 
was  then  decided  to  make  the  value  of  F,  at  any  one  iteration, 
depend  on  the  difference  between ^1T  and  |CSF|  .  An  error  was 
defined  by 

ERROR  =  j  ,n  -  !  csf!  j  (3.1) 

where  the  sum  extends  over  the  entire  PSF  plane.  A  series  of 
different  F  functions  were  tried,  as  defined  by  Table  3.3.  In 
all  of  these,  the  smaller  the  error  becomes,  the  larger  the 
value  of  F  is.  In  Cases  15  to  21,  increasing  the  strength  of 
F  resulted  in  faster  convergence.  In  Cases  22  through  26,  with 
stronger  F's,  the  error  decreased  to  a  certain  point  and  then 
began  to  rise.  A  final  error  of  less  than  20  is  indicative  of 
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a  very  good  estimate  of  the  phase  aberration.  Although  F 
function  XII  worked  well  with  random  number  seed  A,  in  Case  27 
with  random  number  seed  B,  convergence  was  not  obtained. 

Cases  28  to  32  are  with  a  larger  amount  of  aberration. 
Nothing  worked  for  these  cases,  including  a  very  large  number 
of  iterations,  and  a  very  strong  F  function,  which  would  no 
doubt  have  produced  an  instability  in  the  previous  cases. 

Case  4  was  repeated  with  1000  iterations,  which  is  nearing  the 
limit  of  reasonable  computation  times  on  our  11/70,  with  no 
improvement.  A  final  algorithm  was  tried  in  which  every 
time  the  error,  as  defined  by  Equation  (3.1)  began  to  increase, 
a  new  random  number  set  between  0  and  2n  was  added  to  the 
current  estimate  of  W  and  the  algorithm  was  allowed  to 
continue.  After  1000  iterations  of  this,  good  convergence 
was  still  not  obtained. 

The  c;.  i>.  algorithm  is  quite  attractive  from  the  standpoint 
of  implementation  into  Special  Purpose  hardware,  and  the 
overall  simplicity  of  the  method,  but  unless  the  problem  of 
false  convergence  is  solved,  it  will  not  be  useable. 


SI 


r  tor 


Table  3.3  Parameters  for  Convergence  Acceleration 


4.0  THE  DEVANEY  ALGORITHM 


4.1  NOISE-FREE,  POINT  DETECTOR  RESULTS 

The  definition  of  the  Devaney-Gons laves  phase -retrieval  al¬ 
gorithm  is  presented  in  the  proprietary  Appendix  A.  In  this 
chapter  we  merely  present  results.  The  results  are  repeated  in 
the  Appendix  with  additional  discussions  concerning  the  algorithm. 

Figure  4.1,  4.2,  and  4.4  show,  respectively,  contour 
plots  of  the  original  simulated  phase  aberration  across  an 
aperture,  the  aberration  estimated  by  the  program,  and  the 
residual  aberration  (the  initial  as  corrected  by  the  estimate). 
The  values  at  the  contour  lines  are  indicated  by  the  symbols  in 
the  table  to  the  riqht  of  the  graph.  All  values  are  in  terms  of 
waves.  The  aperture  coordinates  X  and  Y  are  in  terms  of  the 
reduced  coordinates.  The  aperture  itself  is  the  best  circle 
which  fits  inside  the  square  array;  since  the  square  array  is 
sampled  on  only  a  9x9  mesh,  the  aperture  is  not  a  perfect  circle 
On  Figure  4.4,  there  is  little  residual  aberration;  less  than  a 
deviation  of  0.1  waves  to  either  side  of  zero,  which  is  the 
minimum  contour  step  value  used  by  the  contouring  program. 

Figures  4.3  and  4.5  show  the  MTF ' s  corresponding  to  the  initial 
aberration  and  the  residual  aberration.  The  latter  MTF  is 
essentially  diffraction  limited;  the  deviation  from  circularity 
near  the  center  is  due  to  the  coarse  grid  on  which  the  MTF  is 
evaluated.  Only  half  of  the  MTF  is  plotted  since  it  is  symmetric 
under  the  transformation  (x,  y) — -y)  .  The  coordinates  x 
and  y  on  the  graph  represent  normalized  spatial  frequencies;  a 
value  of  1  indicates  a  frequency  of  IF*.  Figures  4.6  through 
4. 10  show  contour  plots  of  the  results  for  another  case.  Yet 
another  case  is  presented  in  Figures  4.11  through  4.15.  This 
case  had  larger  aberrations  and  took  longer  to  complete  than  the 
previous  cases. 
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Figure  4.9  Contour  Plot,  of  the  Residual  Phase  Aberration 
(initial  phase  aberration  less  the  estimated  phase  aberration) 
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Figure  4. IS  Contour  Plot  of  MTF 
from  the  Residual  Phase  Aberration 

The  next  case  used  a  larger  and  more  complex  aberration 
and  required  still  more  computer  time  for  proper  performance. 
Figure  4.16  numerically  displays  the  phase  aberrations  over  the 
aperture  in  units  of  radians.  As  can  bo  seen,  the  maximum 
residual  aberration  at  any  point  is  less  than  0.03  waves. 
Figures  4.1?  through  4.19  show  contour  plots  for  this  case;  the 
residual  phase  aberration  and  MTF  arc  not  displayed  since  they 
arc  effectively  diffraction  limited. 
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Phase  Aberration 
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Anothei  case  is  shown  in  Figures  4.20  through  4.22.  Figure 
4.20  is  the  initial  point  spread  function.  The  results  of  apply¬ 
ing  the  wavefront  correction  as  estimated  by  the  phase  retrieval 
program  is  shown  in  Figure  4.21.  This  estimate  was  clearly  not 
an  optimum  estimate.  Nevertheless,  it  is  close  enough  to  the 
actual  aberration  that  the  corrected  point  spread  function  is 
considerably  improved.  The  phase  retrieval  program  was  then  run 
again,  with  this  improved  point  spread  function  as  input. 
Correcting  the  system  this  time  yields  diffraction  limited  perfor¬ 
mance.  This  is  an  illustration  of  an  "iterative"  correction 
scheme,  where  each  correction  results  in  an  improvement  over  the 
previous  state  of  the  optics.  If  the  aberrations  arc  severe 
enought,  it  may  well  not  be  feasible  to  do  all  the  correction  at 
once,  but  rather  to  employ  a  series  of  correction  steps. 


st  Correction 
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EFFECTS  OF  NOISE 


Wo  now  present  results  showing  the  effects  of  detector 
noise.  This  noise  is  white,  Gaussian,  signal- independent  noise 
characterized  by  a  Noise  Equivalent  Flux  Density  (NEED).  This  is 
the  power  from  a  point  source  falling  on  the  entrance  pupil 
which  gives  a  signal  to  noise  ratio  of  unity  in  the  detector 
which  receives  the  maximum  amount  of  power  from  the  dil fraction 
limited  point  spread  function.  As  such  defined,  the  NF.FD  is 
dependent  on  the  area  of  the  entrance  pupil.  Assume  we  are 
sampling  at  the  Nyquist  rate  with  square  detectors  of  area  a,  and 
s.  is  the  star  imago  at.  the  j  ^sample  position.  The  HMS  noise 
powe  r  i s 


o  -  A  '  NEED 

where  A  is  the  aperture  size.  The  total  power  from  one  detector 
is 


where  n  is  a  unit  variance  white  Gaussian  noise  signal.  The  sum 
of  the  power  in  the  noi »o  signal  must  equal  the  incident  power: 


a 


o 


where  1  is  the  power  pet  unit  area  over  the  entrance  pupil.  If 
a  normalized  signal  is  defined 


a 


then 


l 


and 

A 

S  .  =  s .  '  I  A  ♦  A  •  NEFD 
IDO  n 

The  signal  to  noise  ratio  at  the  diffraction  limited  PSF  peak, 

A 

Sp,  is 

!o  * 

<;  /s*  ~  <; 

w  *  NKH)  ‘  p 

\ 

In  our  simulation,  Sp  was  equal  to  69.Z289.,  so  for  a  S/N  ratio 
of  W,  wo  have 


o 

NKFD 


W 


- 1  ‘>  2 
x  1 0  wa 1 1  s/ cm  i‘ 


A  zero  magnitude  star  has  a  power  of  5 
Assuming  a  narrow  band  filter  transmitting  0 .  i.  „  wavelength  to  get 


nearly  monochromatic  radiation,  the  total  power  is  .1.6  x  10 

*) 

watts/cm".  Ttie  star  magnitude  simulated  is  thus,  from  the 
definition  of  magnitude. 


■16 


2 . 4  log . 


i  0 


1:2 '  Kr;rn 

l  <f 


6  x  ’.0 


Figures  4.71  through  4.11  present  the  results  of  a 
simulation  with  noise.  Figure  4.21  shows  contours  of  the  initial 
phase  aberration  across  the  exit  pupil.  The  key  to  the  right 
shows  the  contour  values  in  terms  of  waves.  Figures  4.24  shows 
contours  of  the  MTV  corresponding  to  this  phase  aberration. 

Figure  4.21  is  the  MTK  after  noise  has  been  added .  The  phase 
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Phase  Aberration 


Figure  4.26  Estimated  Phase  in  Noisy  Case 


figure  4.30  Original  PSF  with 


aberration  as  estimated  by  the  program  is  displayed  in 
Figure  4.26.  The  residual  phase,  equal  to  the  initial  less 
the  estimate,  is  shown  in  Figure  4.27.  Note  that  it  is 
significantly  less  aberrated  than  the  initial  phase.  The  MTF 
corresponding  to  the  residual  phase  aberration  is  in  Figure  4.28. 
No  noise  has  been  added  to  this  MTF.  The  point  spread  function 
from  the  initial  aberration  is  shown  in  Figure  4.29.  The  same 
PSF  with  added  noise  is  shown  in  Figure  4.30.  The  PSF  from  the 
corrected  system  is  displayed  in  Figure  4.31.  This  is  definitely 
sharper  than  the  first  PSF.  When  the  same  case  wa3  run  without 
noise,  the  corrected  PSF  was  as  shown  in  Figure  4.32.  Thi  is 
very  nearly  diffraction  limited.  However,  it  must  not  be 
thought  that  the  noise  was  really  responsible  for  a  poore j  phase 
estimate,  when  the  program  was  rerun  in  such  a  manner  as  '.o 
simulate  the  results  of  a  more  complex  version  of  the  algv  ithm, 
a  much  better  estimate  of  the.  phase  aberration  was  obtainc-  . 

The  PSF  from  this  correction  is  displayed  in  Figure  4.33.  The 
version  of  the  algorithm  needed  to  insure  good  results  in  this 
case  requires  much  more  computer  time  to  execute.  For  a 
different  phase  aberration  or  an  algorithm  differing  in  its 
details,  it  could  well  be  that  a  noise-free  simulation  might 
yield  a  poorer  estimate  than  a  simulation  with  noise.  The 
algorithm  is  not  perfect,  but  the  noise  presents  no  additional 
difficulty.  tt:  should  not  be  difficult  to  find  stars  which 
will  provide  at  least  the  signal  to  noise  ratio  tested  here. 

4.3  EFFECTS  OF  NON-POINT  DETECTORS 

In  the  case  of  extended  (non-point)  detectors  the 
measured  point  spread  function  PSF(x,  y)  is  related  to  the  ideal 
point,  spread  function  PSF  (x,  y)  by 

A 

PSF  (x,  y)  *  PSF ( x ,  y)  ®  D ( x , y) , 

whcie  D  (x,  y)  is  the  detector  response  which  we  have  takan 
as  a  square  of  side  I.: 
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Figure  4.32  Corrected  PSF  iror.'  Noise-Free  Case 
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A  detector  size  of  100  units,  as  used  in  the  simulations,  has 
the  first  zero  due  to  the  sine  function  at  a  spatial  frequency 
of  0.01  cyclos/unit.  The  optical  system  with  which  we  have 
simulated  here  has  a  diffraction-limited  cut-off  frequency  as 
determined  by  the  Fi  and  the  wavelength  of  O.O'i  cyclos/unit 
which  is  severely  limited  by  the  detector  size. 

Despite  these  detector  size  limitations  on  the  received 
signal,  the  phase  retrieval  algorithm  shows  a  remarkable  ability 
to  correctly  perform  estimates,  even  in  the  presence  of  noise. 

The  results  obtained  in  one  simulation  are  presented  in 
Figures  4.34  through  4.38.  Figure  34  shows  the  intensity  profile 
of  a  star  object  of  the  same  brightness  relative  to  the  detector 
noise  as  presented  in  Section  4.4  imaged  by  an  optical  system 
jxissessing  a  wave  aberration  function  shown  in  Figure  4..3.r>.  The 
image  of  the  star  was  assumed  to  move  across  the  local  plane 
detector  array  at  a  sufficiently  small  angle  to  the  array  axis 
to  allow  the  convolution  of  the  star  image  with  the  array  trans¬ 
fer  function  to  be  sampled  at  the  Nyqurst  rate,  F*  “2 ,  which  is 
10  units  in  the  simulation.  This  convolution  is  shown  in  Fiqurc 
4.36  whore  it  is  assumed  that  the  detector  ar  ay  is  composed  of 
square,  100  unit  x  100  unit  detector  elements.  The  signal  has 
been  multiplied  by  a  factor  of  ten  relative  to  Figure  4.38,  since 
otherwise  it  would  be  difficult  to  see.  Assuming  the  same  NEFD 
as  vised  in  Section  4.?  the  noise  corrupted,  detected  signal  which 
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Contour  values  fin  waves 


5.0  THE  TWO- IMAGE  PLANE  ALGORITHM 


5.1  INTRODUCTION 

The  two-image  plane  algorithm  was  devised  in  order  to 
perform  wavefront  aberration  estimation  from  an  unknown,  extend¬ 
ed  object.  When  imaging  a  delta-function-like  object,  the 

transform  of  the  measured  signal  S{x,y)  is  P(f  ,f  ),  the 

x  y 

OTF  of  the  imaging  system.  An  extended  object  yields  a  transform 

M(f  ,  f  )  ~  0(f  ,  f  )  •  P(f  ,  f  ) ,  where  O  is  the  object 
X  )  X  y  x  y 

transform.  Since  O  is  assumed  unknown,  P  cannot  be  separated 
out  t.o  permit  use  of  the  one-image  algorithm  already  discussed 
in  Chapter  4.  Appendix  D  presents  the  actual  algorithm  used 
and  provides  a  more  complete  discussion  of  results  than  in  this 
Chapter . 


5.2  RESULTS 

All  the  results  presented  here  used  two  image  planes 
separated  by  one  wave  of  defocus,  that  is,  a  paraboloidal 
aberration  which  is  zero  at  the  center  of  the  aperture  and  equal 
to  one  wave  at  the  aperture  edge  is  added  to  the  existing 
aberration  to  simulate  the  result  of  moving  the  imago  plane.  The 
proper  amount  to  shift  the  image  plane  is  presumably  an  amount 
which  causes  a  significant  change  in  the  OTF;  it  the  do focus  is 
small  compared  to  the  aberrations  already  present,  we  would 
expect  poorer  results.  Variations  in  the  amount  of  shift  were 
not  investigated  here.  In  an  imaging  system  using  a  moveable 
image  plane,  the  amovint  of  de focus  can  be  controlled  as 
appropriate  for  the  aberrations  present. 

Results  are  shown  for  only  one  particular  phase  aberration, 
although  it  is  thoroughly  investigated.  The  initial  RMS 
wavefront  error  was  .708  waves.  All  cases  presented  only 
simulated  the  operation  of  the  extremely  lengthy  "complete" 
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algorithm  which  is  virtually  guaranteed  of  finding  a  good  esti¬ 
mate  in  the  absence  of  noise.  Table  5.1  presents  results  for 
varying  amounts  of  noise.  Point  detectors  and  point  objects  were 
used.  The  sampling  is  assumed  to  be  at  the  Nyquist  rate.  The 
RMS  error  of  the  system  after  correction  by  the  estimate  is  shown 
in  the  column  labelled  "RMS" .  As  the  noise  increases*  -the  error 
of  the  phase  retrieval  estimates  also  increases. 


RMS 

Error 

0 . 708 
■ - 


TABLE  5.  1 

Results  of  the  Two  Image  Plane  Simulations 


Noise 

Level 


Detector- 

Size 


Object 

Width 


Case 


Original 


0.0 

NO’ 

Point 

L  _ _ 

Point 

K 

0.016 

0.01 

Point 

Point 

N 

0,056 

0 . 02 

Point 

Point 

P 

0.081 

o.o-; 

Point 

Point 

0 

Tabic  5.2  shows  simulat '  with  non-point  detectors  (but 
still  point  objects).  As  j  c  one  imago  simulation"'  the 
detectors  arc  assumed  t  '  c  s<;  res  100  units  on  a  side,  with 
a  >  F*  product  of  20  uni  -  .  error  in  the  phase  estimate  has 

increased  by  a  facto*  >.  for  case  U  versus  case  N,  both 

having  a  noise  lev*.  .  .  Case  2*  with  1/8  as  much  noise  as 

case  N  still  has  a  larger  error.  Figure  5.1  shows  MTF's  that 
reveal  how  much  information  is  lost  by  use  of  100  unit  detectors. 
Figure  5.1(a)  is  the  MTF  from  the  initial  aberration  in  case  N. 
Figure  5.1(b)  has  had  the  noise  from  case  N  added.  Figure  5.1(c) 
presents  the  transform  of  the  detected  signal,  which  is  the  MTF 
from  Figure  5.1(a)  multiplied  by  the  transform  of  the  detector 
response.  Note  that  little  information  is  present  0)  ept  at  the 
very  lowc  ..  frequencies.  Figure  5.1(d)  shows  the  result  of  adding 
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the  noise  in  ease  Z  (which  is  much  lens  than  in  case  N) .  These 
MTF's  display  only  one  quadrant. 


TABLE  5.2 

Results  of  the  Two-Image  Plane  Simulations 


RMS 

Noise 
Level  ... 

Detector 

Site 

Object 

Width 

Case 

0.708 

Original 

0.087 

0.01 

100 

Point 

U 

0.070 

0.005 

100 

Point 

V 

0.045 

0.0025 

100 

Point 

Y 

0.022 

0.00125 

100 

Point 

Z 

0.01ft 

0.0)7 

0.338 

0.025 

0.014 

0.000 

0.000 

0.000 

0.000 

0.074 

0.054 

0.037 

0.014 

0.007 

0,009 

0.000 

0.000 

0.000 

0.047 

0.0)7 

0.070 

O.Q77 

0,023 

0.007 

0.020 

0.000 

0.000 

0.105 

0.083 

0.052 

0.048 

0.08b 

0.024 

0.060 

0.018 

0.000 

0.057 

0.060 

0.067 

0.108 

0.053 

0.077 

0.059 

0.020 

0.014 

0.1)6 

0.020 

0.027 

0.142 

0.0)6 

0.071 

0.045 

0.0)4 

0.013 

0.094 

0.046 

0.074 

0.165 

0.027 

0.040 

0.029 

0.082 

0.038 

0.242 

0.024 

0.081 

0.043 

0.028 

0.04) 

0.055 

0.058 

0.035 

1.000 

0.170 

0.102 

0.093 

0.108 

0.1  11 

0.0)2 

0.016 

0.027 

Figure  5.1(a)  Case  N,  Original  MTF 


0.031 

0.059 

0.026 

0.041 

0.049 

0.068 

0.084 

0.038 

0.030 

0.061 

0.093 

0.039 

0.052 

0.025 

0.021 

0.055 

0.092 

0.022 

0.080 

Ci.  05  3 

0.048 

0.114 

0.056 

0,058 

0.089 

0.010 

0.036 

0.070 

0.040 

0.041 

0.015 

0.070 

0,076 

0.055 

0.044 

0.018 

0.077 

0.088 

0.117 

0.083 

0.068 

0.064 

0.032 

0.038 

0.058 

0.143 

0.072 

0.0)4 

0.151 

0.042 

0,005 

0.031 

0.011 

0,072 

0.198 

0.060 

0.063 

0.190 

0.009 

0.045 

0.008 

0.076 

0,038 

0.263 

0.042 

0.061 

0.078 

0.019 

0,060 

0.066 

0.022 

0.048 

1.000 

0.167 

0.095 

0.069 

0.110 

0.135 

0.094 

0.036 

0.033 

Figure  5.1(b)  Case  N,  Noisy  MTF 
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0.001  0.001  0.000  0.000  0.000  0,000  0.000  0.000  0.000 

0.002  0.001  0.000  0.000  0.000  0.000  0.000  0.000  0,000 

0.004  0.002  0.001  0.001  0.000  0.000  0.000  0.000  0.000 

0.008  0.003  0.000  0.001  o.oot  o.ooo  0.000  0.000  0.000 

0.005  0.003  0.001  0.002  0.000  0.001  o.ooo  o.ooo  0.000 

0.023  0.002  0.000  0.004  0.001  0.001  0.001  0,000  0.000 

0.009  0.003  0.001  0.003  0.000  O.OOO  0.000  0,000  0.000 

0.137  0.008  0.004  0.004  0,001  0.00?  0,003  0.00}  0.001 

1.000  0.096  0.010  0.015  0.010  0.0011  0.003  0.000  0.002 

Figure  5.1(c)  Case  ,  Transform  of  Detected  Signal  of  Point  Object. 

0.003  0.004  0.003  0.006  0.008  0.0OH  0.010  0.005  0.004 

0.005  0.006  0.000  0.004  0.004  0.00?  0.007  0,011  0.003 

0.000  0.001  0.003  0.008  0.004  0.007  0.011  0.001  0.004 

0.003  0.003  0.00?  0.005  0.007  0.006  0.003  0.004  0.002 

0.006  0.007  0.005  0.008  0.004  0.003  0.006  0.002  0*007 

0.023  0.005  0.001  0.006  0.006  0.004  0.003  0.005  0.0C9 

0.00b  0.004  0.005  0.005  0.003  0.001  0.003  0.003  0.002 

0.136  0.010  0.008  0,004  0.00?  0,006  0.002  0.006  0.004 

1.000  0,095  0.013  0.019  0.010  0.010  0.008  0.004  0.004 

Figure  5.1(d)  Case  7.,  Transform  of 
Detected  Signal  of  Point  Object  Plus  Noise 


Cases  AA,  Ml,  and  CC  in  Table  5.3  have  the  same  noise  and 
detector  size  as  Case  7.,  but  now  an  extended  object  is  used.  The 
object  is  taken  to  have  an  intensity  in  the  image  plane  given 
by  a  Gaussian 

IOBJ(*-1>  *  CXP  *  *2>''DOBjl- 


TABLE  5.3 


Results  of  the  Two-Image  Plane  Simulations 
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t>0B7  is  listed  in  the  table  under  "Object  Width"  and  is  in  the 
same  units  as  the  detector  size.  The  detected  signal  is  given 
by  the  convolution  of  the  system  1>SF  with  the  detector  response 
and  the  object  distribution: 

S(x,y)  =  PSF(s,y)  •  D(x,y)  •  I0BJ(x#y) 

The  transform  of  IQBj  is 

WV  V  -  '*>>  [-  ”2r>0BJ,fx  +  fy»] 

The  transform  of  the  detected  signal  is 

s(fx'  V  “  fy>  D(£X-  V  •  Wfx-  V 

As  the  object  size  increases  (in  case  AA,  BB,  and  CC  in  Table 
5.3),  the  error  on  the  phase  retrieval  estimate  increases.  The 
additional  attenuation  of  the  high  frequencies  due  to  the  object 
size  is  shown  in  Figure  5.2,  which  represents  the  detected  signal 
transform  with  and  without  noise,  from  case  CC.  This  may  be 
compared  to  case  7.  in  Figure  5.1(c),  in  which  the  only  difference 
is  the  use  of  a  non-point  object  in  case  CC.  More  high  frequency 
information  is  lost  in  case  CC.  Case  DD,  from  Table  5.3  has 
twice  as  much  noise  as  case  CC,  and  may  be  compared  with  case  Y, 
which  has  the  same  amount  of  noise  but  a  point  object.  Case  DD 
yielded  an  estimate  having  twice  the  residual  error  of  case  Y. 
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0.000 
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O.oon 
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0.000 

0.000 
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0.000 
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0.000 

0.000 

0.000 

0,040 

0.001 

o.oon 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

1  .000 

0.028 

0.000 

0.000 

0,000 

0.000 

0.000 

0.000 
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Figure  5.2(a)  Case  CC, 

Transform  of  Detected  Signal  of  Extended  Object 
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0.001  0.001  0*003  0.006  n.GOH  0.008  0.010  O.OOS  0.004 

0.003  0.00k  0.009  0  *  004  0.004  0.00?  0.00?  0.011  0.003 

0.004  0.002  0.003  0.008  0.004  0.007  O.OU  0.001  0.004 

O.OOS  0,006  0.00?  0.004  0.000  0.006  0.003  0.004  0.00? 

0.004  0.004  0.006  0.006  0.004  0,004  O.OOS  0.002  0.0C7 

0.004  0.000  0.001  O.OOS  0.006  0.004  0.004  0.005  0.004 

rt.012  0.003  0.005  0.00?  0.003  0.001  0.003  0.003  0.002 

0.041  0,003  0.00.3  O.OOS  0.C03  0.006  0.001  0.006  0.004 

1.000  0,077  0.004  0.003  0.001  0^002  0^009  0.004  0.005 

Figure  5.2(b)  Case  CC, 

Transform  of  Detected  Signal  of  Extended  Object  Plus  Noise 

In  summary,  the  two-image  plane  technique  certainly  works; 
it  does  suffer  from  the  same  problem  as  the  one-image  plane 
algorithm,  that  of  very  long  execution  times  to  insure  the  best 
estimate.  Finally,  it  may  be  worth  noting  that  although  we  have 
used  the  addition  of  defocus  to  the  system,  since  it  is  easy  to 
add  a  known  amount,  any  sort  of  other  aberration  could  also  work, 
as  long  as  it  was  known.  These  aberrations  could  be  induced  by 
varying  the  active  optics,  or  by  the  introduction  of  aberrating 

elements  into  the  optical  system.  Also  possible  are  algorithms 
utilizing  three  or  more  different  focus  positions  which  may  be 
useful  in  the  presence  of  a  large  amount  of  noise,  or  for  very 
large  objects. 


6.0  CONCLUSIONS  AND  RECOMMEND AT IONS 


EIKONIX  has  successfully  demonstrated  phase  retrieval  of 
a  complex  amplitude  from  its  modulus  under  a  wide  variety  of 
conditions.  Our  methods  are  directly  applicable  to  the 
problem  of  actively  controlling  an  optical  system  and  possess 
definite  advantages  over  other  methods  such  as  direct  inter¬ 
ferometric  measurements  of  the  wavefront  aberration  which 
requires  integration  of  elaborate  and  costly  electromechanical 
devices  into  the  system.  In  particular,  our  algorithm  works 
well  for  a  wide  variety  of  wavefront  aberrations,  reference 
objects,  detector  sizes  and  noise  levels. 

The  only  information  needed  to  perform  phase  retrieval  is 
the  detected  signal  of  a  bright  "point"  object  sampled  at  the 
Nyquist  rate.  Although  the  planned  detector  sires  would  seem 
to  preclude  Nyquist  sampling,  use  of  image  motion  across  the 
focal  plane  array,  with  suitably  spaced  time  samples,  can 
provide  the  necessary  spatial  resolution.  Despite  the  fact  that 
much  of  the  higher  frequency  information  in  the  OTF  is  lost  due 
to  the  detector  size,  our  simulations  were  successful  with  the 
planned  HALO  detector  sizes. 

Phase  retrieval  simulations  have  not  been  performed  for 
polychromatic  radiation  but  OTF  calculations  show  that  for  a 
bandpass  filter  transmitting  from  0.8  to  1.2  of  t.ho  central 
wavelength,  the  polychromatic  OTF  differs  only  slightly,  and  at 
the  higher  frequencies,  from  the  monochromatic  OTF.  This 
difference  would  be  obscured  by  noise  and  the  large  detector  size. 
Thus,  we  do  not  feel  that  there  would  be  any  difficulty  in  per¬ 
forming  phase  retrieval  on  broadband  radiation  provided  some 
sort  of  filter  is  provided. 

Simulations  have  been  run  using  white  noise  on  the 
detected  signal  as  specified  by  the  NKFD  of  the  focal  plane 
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detector  arrays.  Phase  retrieval  has  been  shown  to  be 
possible  even  when  visually  the  signal  from  a  point  source  is 
apparently  totally  obscurred  by  the  noise.  The  noise 
simulations  also  assumed  a  very  narrow-band  spectral  filter 
which  reduced  the  energy  from  the  point  source  significantly. 

It  should  not  be  difficult  to  find  stars  which  yield  a  sufficient 
signal- to-noise  ratio. 


It  may  be  necessary  to  estimate  aberrations  in  regions  of 
the  field  of  view  which  are  totally  below  the  horizon,  and  con¬ 
sequently  will  probably  lack  bright  point  objects.  In  this 
case,  it  is  still  possible  to  perform  phase  retrieval  using  a 
bright  extended  (non-point)  object  as  a  target;  however,  two 
images  must  bo  obtained  at  different  focus  positions.  This 
procedure  has  been  successfully  simulated. 

The  chief  difficulty  with  the  method  lies  in  executing 
the  algorithm  in  an  acceptable  length  of  time.  Current 
simulations  take  about  thirty  minutes  on  a  POP  11/70.  Special 
purpose  hardware  could  reduce  this  by  at  least  a  factor  of  100 
but  the  computation  time  increases  extremely  rapidly  as  the 
magnitude  and  the  number  of  degrees  of  freedom  of  the  wavefront 
aberration  increases.  It  will  be  necessary  to  reduce  the 
running  time  before  realistic  simulations  of  the  HALO  system 
could  be  performed. 

Further  study  of  the  Gerchberg-Saxton  algorithm  is  needed 
to  understand  its  convergence  problems  in  two  dimensions.  * t 
works  well  in  one  dimension  and  has  the  potential  to  be  much 
faster  than  the  parameter  search  method. 

Another  question  to  be  answered  is  the  range  of  magnitudes 
of  the  aberrations  such  that  phase  retrieval  can  be  performed  in 
a  specified  period  of  time.  This  will  impact  on  how  accurate 
the  initial  alignment  procedure  must  be.  This  question  must 
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also  bo  coordinated  with  a  more  realistic  model  of  the 
aberrations  expected  from  HALO,  both  at  the  point  of  takeover 
of  control  after  the  initial  alignment,  and  the  drifts  in  mirror 
positions  as  a  function  of  time.  Also  to  be  considered  are  the 
effects  of  short  term  fluctuations  which  could  cause  shifts  in 
the  aberrations  between  measurements  done  at  slightly  different 
times,  which  otherwise  would  have  been  expected  to  be  measuring 
exactly  the  same  aberration  functions. 

The  problem  of  precisely  which  stars  that  arc  bright 
enough  to  be  used  as  print  ob]ects  and  which  will  be  able  to  be 
acquired  is  not  yet  answered.  Although  there  are  many  stars 
which  would  be  suitable,  there  is  the  additional  constraint  that 
nearly  simultaneous  measurements  be  madr  in  perhaps  six  different 
isoplanatic  regions. 

Simulation  of  extended  (non-point)  objects  has  been 
limited  to  Gaussian  intensity  distributions.  Other  kinds  of 
objects,  such  as  random  scenes,  may  be  more  appropriate. 

The  problem  of  “deconvolution",  that  is,  obtaining  mirror 
positions  from  a  knowledge  of  the  wavefront  error  in  several 
different  regions  of  the  field  of  view  is  complicated  by  the 
two-fold  ambiguity  in  the  phase  retrieval  estimates  inherent  in 
any  problem  with  a  symmetrical  pupil  function.  Although  many  of 
the  field  points  will  not  have  a  nearly  synuuetr ical  pupil,  the 
magnitude  of  this  problem  must  be  determined  in  conjunction  with 
the  designers  of  the  HALO  optical  system.  Simulations  so  far 
have  only  treated  the  case  of  on-axis  imaging  of  a  symmetric 
optical  system,  Eventually  off-axis  cases  must  be  considered, 
and  possibly  non-symmetr ic  systems.  A  realistic  treatment  will 
involve  more  detailed  knowledge  of  the  optical  system. 
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